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Abstract — This paper proposes a new framework for the 
optimization of excitation inputs for system identification. The 
optimization problem considered is to maximize a reduced Fisher 
information matrix in any of the classical D-, E-, or A-optimal 
senses. In contrast to the majority of published work on this 
topic, we consider the problem in the time domain and subject to 
constraints on the amplitude of the input signal. This optimization 
problem is nonconvex. The main result of the paper is a convex 
relaxation that gives an upper bound accurate to within 2/tt 
of the true maximum. A randomized algorithm is presented for 
finding a feasible solution which, in a certain sense is expected to 
be at least 2/ty as informative as the globally optimal input signal. 
In the case of a single constraint on input power, the proposed 
approach recovers the true global optimum exactly. Extensions 
to situations with both power and amplitude constraints on 
both inputs and outputs are given. A simple simulation example 
illustrates the technique. 

I. Introduction 

System identification is the process of computing a com- 
pact mathematical model of a real-world system based on 
experimental input-output data. The quality of the model so 
identified can depend a great deal on the choice of excitation 
input. In many practical applications it is natural to seek 
to extract as much relevant information from the system 
as possible in minimal time, subject to certain experimental 
constraints. 

Over several decades, a large body of literature has devel- 
oped on the topic of optimal input design (see, e.g., (TJ, 0, 
[3 1, [4], [5 1, |6| and references therein). Essentially the same 
problem is studied in the communications literature for finding 
test signals for channel estimation (see, e.g., [7], |8|, |9| and 
many others). Most channel estimation systems assume quite 
simple dynamic models - either a static input-output map or 
an FIR filter - however recently more dynamic models have 
been considered [ 10 1 . 

The bulk of input design methods for dynamic systems are 
based on the recognition that, for a linear system, the Fisher 
information matrix is an affine function of the input power 
spectrum. Imposing power constraints via Parseval's theorem, 
and optimizing over an affine parametrization of the input 
spectrum, the optimization can be posed as a semidefinite 
program, for which efficient computational tools are readily 
available (3). 

In the frequency domain it is natural to consider signal 
power consttaints. However, in many practical cases, the real 
constraint is the amplitude of the excitation input, not its 
power. In industrial processes, amplitude constraints are com- 
mon, as evidenced by the success of model model predictive 



control. Furthermore, in the emerging area of biomedical 
system identification safety limits are often given as amplitude 
constraints (see, e.g. ifTTI . Ifl2ll ): in communication channel 
estimation, often a binary signal is desired. The relationship 
between signal phases and peak amplitude is highly non- 
linear and non-smooth, making optimization under amplitude 
constraints computationally challenging in the frequency do- 
main. In previous work we have applied a Pdlya-like algorithm 
to this problem, solving first a convex optimization with 
power constraints followed by a sequence of smooth nonlinear 
optimizations 1 1 3 1 . 

In the time domain, amplitude constraints appear naturally, 
and one can study linear time-varying systems and time- 
varying constraints, as may be appropriate to estimate intrinsic 
parameters of a nonlinear system via small deviations about a 
changing operating point. However, the resulting optimization 
problem is highly nonconvex: even for a system with one 
parameter to identify, the optimization for an input signal of 
length n will have 2" local optima. Recent work has suggested 
using a frequency-domain power-consttained optimization as 
an initial guess for a local BMI optimization algorithm |14|. 
Others have suggested a minimizing one-step-ahead parameter 
variance [15], or optimizing the transition probabilities of a 
markov process for the input fl6l . In this paper, we show 
that techniques from semidefinite relaxation of non-convex 
quadratic programs can be extended to the time-domain ex- 
periment design problem to find approximate solutions with a 
high degree of efficiency. 

A. Semidefinite Relaxations of Nonconvex Quadratic Pro- 
grams 

As shall be seen, the problem of maximizing an information 
matrix in the time domain has a sttucture similar to certain 
nonconvex quadratic programming problems. Such problems 
are in general NP-Hard to solve exactly, however it has been 
found that specific semidefinite relaxations can give upper 
bounds on the objective and lead to efficient randomized 
algorithms to find feasible suboptimal solutions with a high 
degree of accuracy. The breakthrough result of Goemans and 
Williamson IfTTI for finding the maximum cut of a graph to 
within approximately 0.87 of its true optimum led to many 
more applications in combinatorial optimization, signal pro- 
cessing, operator theory, and systems analysis [18|, [19], |20|. 
An essential tool in this paper will be Nesterov's extension 
to maximization of a positive-definite quadratic form over 
a hypercube with an accuracy ratio of 2/7T ED . Ifl8l . This 
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family of methods can be variously interpreted as either a 
simple relaxation of the feasible-set, the dual of the Lagrangian 
dual, or optimization of the covariance of a random variable 

cdd, ma. 



B. Paper Structure 

The Structure of the paper is as follows: in Section [II] we 
introduce the problem statement mathematically; in Section[lIl| 
we give a convex relaxation of the input design problem with 
input amplitude constraints, and the main theoretical results 



of the paper; in Section IV we extend the solution to more 



general contraint types; in Section [V] we give some illustrative 
examples; Section VI contains some brief conclusions and 
future directions. 



C. Notation 

the cone of symmetric n x n positive-semidefinite 
matrices. For symmetric matrices, X > Y means X— Y £ <!?+. 
The function sgn(-) : E" — > { — 1,0,1}™ computes a vector, 
each element of which is the sign of the corresponding 
element of the argument vector. The symbol E denotes the 
expectation operator. We make the following definition: a 
function v : S£ — > K is denoted nonnegative-concave if 
v{X) >0VIe S+ and v(aX +(l-a)F) > av(X) + (l- 
a )v(Y) yx,Y e S+, a e [0,1]. 

II. Optimal Experiment Design 

In a statistical experiment design, the amount of information 
about parameters 8 contained in the observations y from an 
experiment is measured by the Fisher information matrix 1(6), 
which depends on the experimental conditions. The Fisher 
information matrix is defined as 



Is :=E 



d\ogp(y\8) d\ogp(y\8)' 
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where p(y\8), considered as a function of 8 with fixed observa- 
tions y, is the likelihood function, and the derivitives are taken 
at the true value of 8. The inverse 7(0) _1 is a lower bound 
on the achievable covariance matrix of an unbiased estimator 

nza. 

We consider dynamic system estimation problems, where 
from observations of finite-length sequences u(t) and y(t) one 
must estimate a system: 

y(t)=G e u(t)+H e e(t) 

where Qg and Hg are unknown linear maps parametrized by 
8 and e(t) is a Gaussian white noise sequence. Note that this 
framework naturally allows multi-input multi-output systems 
via stacking inputs and outputs, as well as time-varying linear 
systems. 

For simplicity we address here the particular case where 
Qe = G{q) and He — V.(q) are single-input single-output 
finite-dimensional LTI systems given as rational functions of 



the shift operator q. The more general cases are a straightfor- 
ward extension of our method. With this structure, the log- 
likelihood function is given by: 



log p{y\8) 
where 



■ log(27r) 



log(cr e ) 



1 

— y 

e t=i 



e(tf 



e(t) :=H e {q)- 1 [y{t)-G e {q)u{t)}. 

Assuming an open-loop experiment, zero correlation between 
u(t) and e(t), and independently parametrized system and 
disturbance model - i.e. 8 = [8q 6h]' - then the information 
matrix can be decomposed as 



h{u) 0' 
* 



where Ig(u) is the block corresponding to 8q and depends 
on the input, and * is the block corresponding to 8h and 
depends only on the disturbance, and thus cannot be optimized 
by choice of input. Hence optimizing 1(8) is equivalent to 
optimizing the upper-left block: 

Now, suppose the system model has N components, i.e. 
8q € Mr, and consider 

m) :=-'Ho H (q)- ld ^^uit). 

The system Fg is a linear system with one input and N 
outputs, each output representing the sensitivity of one of the 
parameters in 8q to the choice of input. Let Fi,i = 1, 2, ...N 
denote the rows of Fg(q). Note that 



de(t) 
~d8^ 



7i(g)«(t) 



T N (q)u(t) 



Let us define a stacked control vector u := 
[u{l),u(2),...,u(n)]' and define a matrix G W nxn 
representing the action of each Ti on u, i.e. let fi(t) be the 
impulse response of Ti, then Fi is the Toeplitz matrix: 



Fs := 





Mi) 



fi(n) /i(n-l) 






fill) 



(1) 



Then one can represent the elements of the reduced informa- 
tion matrix as 

n 

Ie(u)i,j = X^(«)u(*))(^(«)u(t)) = u'FiFjU 



t=i 



Hence we have 



u'F[F lU 
u'F' N F lU 



u'F[F N u 



u'F' N F N u 



(2) 
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as a compact expression for the reduced information matrix 
in terms of quadratic forms in u. In the general MIMO or 
time-varying cases, an equivalent Ig(u) could be similarly 
constructed by stacking control inputs and computing the time- 
varying equivalent of Fi. 

Note that in most cases Q (q) is nonlinearly parametrized by 
9, and hence the above computations depend on the true value 
of 8. This seems somewhat paradoxical, but in most practical 
cases a reasonable guess for 8 can be made, or multi-stage 
adaptive or robustified optimizations can be performed 0. 

A. Optimality Criteria 

The purpose of input design is to maximize, in some sense, 
the information matrix. In this paper, we consider maximizing 
an objective function of the form 

v(u) = J[I e (u)} 

where J(-) : — > K is any nonnegative-concave function, 
which acts on the reduced information matrix Iq (u) generated 
by u. 

In the experiment design literature the following optimiza- 
tion criteria are common, and all are nonnegative-concave: 

• D-Optimality: Jd[Ig(u)} '■= det[Ie(u)] », 

• E-Optimality: J^[Je(u)] :— mineig[Ig(u)], 
. A-Optimality: J A [I e {u)] := - Tr^u) -1 ]. 

Note that D-Optimality is usually defined as maximizing 
c\ct[Io(u)], however this is not concave and any u achieving 
this clearly also maximizes det[/e(u)]™ which is concave. 
Another possible concave function with equivalent maxima 
would be logdct[/g(u)]. 

B. Constraints 

In the initial part of the paper, we will consider amplitude 
constraints on the input signal. These constraints may be time 
varying: 

\u{t)\ < c(t) Vt = 1,2, ...n, 

for some positive constraint sequence c(t). More general 
constraints including constraints on power and output signals 
will be considered in Section [IV] 

III. Amplitude-Constrained Problem 

The time-domain amplitude-constrained input design opti- 
mization problem can be expressed like so: 



v := max J 

u£R" ,\ui\<Ci 



u'Qn,iU 



u yijvti 



(3) 



where each Qjj g R n >™, each diagonal block Qa > and the 
matrix 

QlA ■ ■ ■ Ql,N 



!N,1 



N 



> 



(4) 



We first note that by making the substitution U = uu' , 
and from the fact that trace is cyclic, for each element of the 



information matrix u'Qi ju = Ti(uu'Qij) = Tr(UQij). In 
this form, the optimization problem ([5]) is equivalent to 



v* = max J 



Tr(£/Q M ) 



Tr(UQ Njl ) ... 
where the feasible set C is defined as 



Tr(UQ ltN )' 
Tr(UQ N , N ), 



C :={UeS+ : U iti < C 2,rank(C/) = 1}. 



(5) 



(6) 



Note that the maximization is now concave in the decision 
variable U but the feasible set is non-convex due to the rank 
constraint: in general, a convex combination of two rank-one 
matrices can have rank two. The constraints that U be positive- 
definite and have diagonal elements less than one are both 
convex. 

Following im . ll2TI . we "relax" this nonconvex optimiza- 
tion problem by dropping the troublesome rank constraint 



Vf/ := max J 

UER 



Tr(C/Qi,i) 



Tr(UQ Ntl ) 
where the relaxed feasible set is 



Tt(J7Qi,jv)' 
Tr(UQ N , N ), 



R:={Ue S+ : £/ M < c}} 



(7) 



(8) 



This optimization problem is the maximization of a concave 
function over the convex cone of semidefinite matrices, subject 
to affine constraints, and hence can be efficiently solved. 

The maximization (|7|l has the same objective function as 
<|3j but a larger feasible set, i.e. C C R. Hence the relaxed 
problem provides an upper bound, i.e. v* < vr. The main 
result of this section is that the upper bound given by the 
convex relaxation is quite tight: 

Theorem 1: The true optimal value v* and the optimal value 
of the relaxed problem vr satisfy the following inequalities: 



-VR <V*< Vr. 



(9) 



To prove this theorem will make use of the following 
theorem of Nesterov on the tightness of SDP upper bounds 
on nonconvex quadratic program: 

Theorem 2: [21 1 Let 

v Qp{Q) = max x'Qx (10) 

x£]R n , \xi I <Ci 

vsdp{Q) = maxTr(QX) (11) 

X £ Ft 



then for any Q € S£ 
2 



□ 

Note that 



vsdp{Q) < vqp{Q) < v S dp(Q). 



vqp(Q) = maxTi(QX). 

-AGO 
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An essential step in the proof of Theorem [2] is the following 
statement which we will also make use of: if x is a Gaussian 
random variable and X is its covariance matrix, then 

2 2 
E(sgnxsgna;') = — arcsin(X) > —X (12) 

7T 7T 

where arcsin(-) of a matrix denotes taking the arcsin elemen- 
twise. 

Proof of Theorem [7} 

Since the relaxed feasible set is strictly larger than the true 
feasible set it is guaranteed that v* < Vr. Therefore we only 
need to prove that 

2 

-VR < V*. 
7T 

Since J is a concave function over the symmetric matrices, 
there exists a representation: 

J{U) = mm[h(i) + Tr(H(i)U)] 

i 

where i varies over a possibly infinite set. Furthermore, since 
J(U) > for all U € S+, then it follows that H{i) > 
and h(i) > for all i, otherwise there would exist a U G S„ 
making h(i) + Tr(H(i)U) negative for some i, and hence 
making J(U) negative. 
Let 

vr; = maxmintfefi) + Tr( H(i)U)] 

uec i 

And let Uc and ic be a matrix and function index for which 
this maximum is achieved. 
Similarly, let 

vr = maxmin[/i(i) + Tr(H(i)U)] 

and let Ur and ir be a matrix and function index for which 
this maximum is achieved. 
By Theorem [2] fixing i = ic, 

2 max[Tr( H(i c )U)} < max[Tr( H(i c )U)} (13) 

7T UGR UEC 

From h(-) > and < f < 1 it follows that f h(i) < h(i) Vi, 
and from this and ( p~3] > we have: 

-rnax[/i(z c ) + Tr( J ff( ic )C/)] 
< max.[h(i c ) + Tr(H(i c )U)\ = v c . (14) 
Now, clearly for all U 

mm[h{i) + Tr(H(i)U)] < h(i c ) + Tr(H(i c )U), 

i 

hence 

Vr = maxmin[/i(i) + Tr(H(i)U)] 

< max[h{ic) + Tr(H{i c )U)] (15) 

U&B, 

and so from (jT3J and ( |T4] > we have 
2 

-Vr < vo- 
ir 

This completes the proof of the Theorem. 



A. Finding a Feasible Solution 

The solution of the convex relaxation (|7]i is an n x n matrix. 
To find an identification input, we need to somehow extract 
a vector of length n. The following randomized procedure 
is common in semidefinite relaxation and has proven to be 
effective in practice |fl9l . 

Compute a matrix D > such that U = D'D Sample a 
vector £ € IR™ with each element an independent normally 
distributed random variable. Compute a candidate solution 



diag(c) sgn(D'£), 



(16) 



where diag(c) is a square matrix with the elements of the con- 
straint sequence c on the main diagonal, and zeros elsewhere. 

Let Ur be the solution of the relaxed optimization problem 
(FT), and Ib(Ur) be the reduced information matrix with U = 
Ur, i.e. 



Ie(U R ) = 



Tr(UQ 1A ) 



Tt(UQn,i) 



Tt(UQ 1iN ) 



Tt(UQ 



N.N J 



then we can state the following theorem: 



Theorem 3: The expectation of the reduced information 



matrix generated by (16 1 satisfies the following bound 



E(T(u)) > -I{Ur). 



(17) 



Essentially, one can say that the randomized strategy is 
expected to give inputs at least 2/tt as informative as the 
solution of the relaxed problem, in terms of the reduced 
information matrix. 

Proof of Theorem ^ 

The proof will use on the following lemma: 
Lemma 1: Let 



A 



Ai. 



JV 



.4 



AM 



A 



N,N 



(18) 



where each block Aij is square of dimension n. Suppose 
A > 0, then A Tl - > where 



A 



Tr 



Tr(A M ) 



Tr(Ai !N y 
Tr{A NtN ) 



(19) 



□ 

Proof: Define vectors ej(x) € K ra for i — 1,2, ..JV to 
have x at the i th element and zeros at every other element, 
e.g. e\(x) = [x • • ■ 0]', consider z £ Mr and consider 
efe(z) to be the vector [ek(zi)' e^z-i)' ■ ■ ■ e\-[z]<i)'\. Then we 
have 

e k (z)'Ae k (z) = z'A k z 
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where Ak is formed by selecting out the (fc, k) element of 
each block Aij, denoted Aij(k,k): 



Ak := 



A x ,i{k,k) 



A NA (k,k) 



AiM(k, k) 



A NtN (k,k) 



Since A > 0, ek(z)' Aek(z) > for any z, so A^ > for 
each k. Now 

= Ai + A 2 + ...An 
which is clearly positive semidefinite since each Ak > 0. ■ 

To prove Theorem [3] we must show that 



E (J(fi)) > -I{U R ) 



Now, 



E (!(«)) - -J(tf fl ) 



'■&(Mi,i) 
Tr(Mjv,i) 



Tr(M 1>JV )' 



(20) 



where 



- Qi,j ( E(flfi') - 2 U R 

7T 



Now, since u = sgn(£) with £ a Gaussian random variable, it 
follows from ( fT2| that 

E (M') - -U R = -(axcsin(J7 fi ) - U R ) > 0. (21) 

7T 7T 

Therefore one can define E > such that 



££" = E(mt') - -C/fl 

and by cyclic trace an equivalent formulation of E (-f(u)) 
~I(U R ) is given by (|20j) with 



M, 



E Qi jE. 



With this substitution, and the fact that the block matrix Q 
is positive semidefinite, it is clear that the block matrix M is 
semidefinite. Now, from Lemma [T] and p0| > it follows that 

E (/(«)) > 1~I(Ur). 
This completes the proof of the Theorem. 



IV. Relaxations for Problems with Power and 
Output Constraints 

In this section we extend the above relaxation approach to 
a broader variety of signal constraints. 



A. Constraints on Input Power 

First we examine a single constraint on the input power: 

IMI2 < Pu 

and show in this special case that the relaxation approach finds 
a global optimum. 

The power-constrained input design problem is to find 



max J 

ueR",|Mli<p„ 



u'Qn,iu 



u'Qn.nu 



(22) 

Using again the substitution U = uv! we define the true 
constraint set, and the relaxed constraint set dropping the rank 
constraint: 

C P = {U e S+ :TrU < p„,rank(ET) = 1}, 
R P = {U eS+ :TrU < Pu }, 



and then ( 22 1 is equivalent to: 

/TT>(tfQ M ) 
v* = max J '. 

\[Tr{UQ Ntl ) 
and the relaxed problem is 

"Tr([/Q M ) 

vr 



max J 

ueiip 



Tr(UQ NA 



Tt(UQi,n) 
Tr(UQ N , N ) 

■ Tr(UQ N , N ) 



(23) 



(24) 

For this special case it happens that the relaxation actually 
attains the optimal value. This theorem is analogous to results 
on the "hidden convexity" of trust-region optimization prob- 
lems (see, e.g., 1123] and many others). 

Theorem 4: Let U* be a solution of the convex optimization 



(24 1. Take u to be the eigenvector corresponding to any 
nonzero eigenvalue of U*, e.g. the largest eigenvalue. Then 
u achieves the global optimum of the nonconvex optimization 
221. 



Proof of Theorem ^ 

Consider the "true" optimization over C: 

U* = argmaxminf/i(i) + TifHUW)] 

uec i 

Fix ic to be the index at which this optimum exists, then we 
can consider equivalently 

maxTrfi? Uc)U) — max u'H(in)u 
uec um n ,u'u<i 

Since H(ic) > we can take the eigendecomposition of 
H(ic) = VAV' with VV = I is orthogonal and consider 
the change of variables z = V'u. Then we have the following 
equivalent optimization problem: 

max z'Az = max zfAu. 

z£W ,z'VV z<l 2gl",z'z<l 

The maximum value of this optimization is the largest diagonal 
element of A, i.e. the largest eigenvalue of H(ic)- Now, 
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consider the same optimization over the relaxed feasible set 
under the change of variables Z = V'UV: 



maxTr(H(i c )VZV')] 



maxTr(V H(i c )VZ)} 

ZeR 

maxTr(AZ)l 
zeR 



max 

zeR 



(25) 



It is clear that optimizing solution Z is a diagonal matrix 
with rank at most k, where k is the multiplicity of the largest 
eigenvalue of A (and hence H(ic)), and this eigenvalue is the 
maximum value of the optimization p5j ). It follows that if u is 
any eigenvector H(ic) corresponding to this eigenvalue, then 

u'H{i c )u = max Tr(H(i c )U) 

U <ZLC p 

= max TT(H(i c )U) 
ueRp 

and any solution of the relaxed problem 

U* = argmaxminf/i(i) + Tr(H(i)U)} 

UeR i 

is a matrix of rank at most k. Take an orthogonal eigendecom- 
position of U and any eigenvector corresponding to a nonzero 
eigenvalue is a solution of the optimization problem ( |22| . This 
completes the proof of the Theorem 



B. General Constraints on Input and Output Signals 

In this section we consider optimizing the information ma- 
trix subject to more general power and asymmetric amplitude 
constraints on both input and output. It will be shown that 
similar relaxation methods can be directly applied. 

Suppose we have an input and output power constraints: 

IMI2 < Pu, \\y\\i < p v 

as well as possibly asymmetric and time-varying constraints 
on the input and output: 

^min,i — — ^niax.ij Vmin.i — Ui — 2/max,j, 

for all i = 1, 2, n. 

Let Gi be the i th row of G, the Toeplitz matrix representing 
the system u — > y, generated similarly to [T] with the impulse 
response of Q. Then with decision variables U € S£ and 
u 6 M. N the above constraints can be put in matrix form: 

Tr((7) < Pu , 
Tr(UG'G) < p y , 



r ^t{UGjGi) GiUymax,i Gi1iy m in,i ^ ?/max,i2/min,Z3 

u = uu'. 

where the amplitude constraints are enforced for all i = 
1,2, ...,n. In this case, the constraint U = uu' is nonconvex, 
however it can be relaxed to the convex constraint U > uu', 
which can be represented in Schur complement form 



With this relaxation, we again have a convex constraint set in 
the variables U and u. 

Note that the output constraints here are for the noise- 
free response of the nominal model. If output constraints are 
strict, a conservative bound may be appropriate to account for 
disturbances and model uncertainty. 

Having found a feasible solution of the relaxed problem, a 
candidate input can be chosen via a randomized strategy: 

u — u + aD't; 

where D'D = U — uu', £ € R™ is sampled from a normal 
distribution, and a > is a scaling parameter chosen as 
large as possible such that u satisfies all the constraints. This 
randomized solution could be further improved using a local 
optimization as in |14|. 

One might ask if uniform approximation accuracy bounds 
can be found as in the case of input amplitude constraints. This 
is unlikely, since it is known that for semidefinite relaxations 
of quadratic maximization subject to multiple ellipsoidal con- 
straints, even in the case where the ellipsoids have a common 
center, the quality of the bound degrades logarithmically in 
the number of constraints ll24l . In the general case described 
in this section there are 2n + 2 constraints, which will be large 
for typical input design problems, and the provable optimality 
bounds are not very optimistic. However, the author has found 
this method to be effective in practice. This will be further 
explored in future publications. 

V. Illustrative Examples 

In this section we show some results on a simple illustrative 
example: 

G(q) = 2 ^ b ~ , H = l, (26) 

with b = 0.1, a-y = — 1.8,a2 = 0.9. The step response of 
this system is shown in Figure [T] We applied the proposed 
algorithm with a constraint that \u(t)\ < 1 for all t and a time 
length of 100 samples. Note that experiment length is quite 
short in comparison to the transient dynamics of the system, so 
asymptotic arguments used in frequency-domain input design 
may not apply. The objective function J(-) := det(-) 1 /™ was 
used. The upper bound on the objective computed via the 
relaxation is 1.82xl0 4 . 

To illustrate that the randomized search for a feasible 



solution (see Section III-Ai is in some sense "intelligent" 



U 
a' 



> 0. 



we generated 50,000 candidate sequences using (16 1, and 
computed the information matrix and the resulting objective 
value. The same process was performed for a purely random 
binary sequence, i.e. sgn(£) where £ is a Gaussian white 
noise sequence. Figure [2] shows a histogram of the ratios of 
objective functions these inputs to the upper bound Vr com- 
puted via the relaxation. It can be seen that almost all signals 
generated using the proposed approach are of high quality 
and one could sample far fewer random candidates, although 
the testing candidates is not computationally expensive. In 
contrast, purely random binary sequences give substantially 
worse results. It is clear that the relaxation-based strategy 
biases the randomization heavily towards good inputs. 



1.8 




Fig. 1. Step response of the example system (26). 
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Optimality ratio 

Fig. 2. Histogram of estimated optimality ratio of random candidate signals 
generated with \16\ as well as purely random binary sequences. A value of 
one would prove global optimality. 



The best value achieved for the objective was 1.54xl0 4 . 
The approximation ratio of the feasible solution to the upper 
bound is 0.85, significantly better than the 2/tt s» 0.64 which 
is guaranteed by Theorem [T] Note that we do not know what 
the true optimal value v* is, only that it is between best 
feasible input found and the upper bound vr. Therefore the 
feasible input we have found may in fact be closer to the 
global optimum than the ratio 0.85 suggests. 

The response of the system to the best input is plotted 
in Figure [3] The usefulness of the input design procedure 
was evaluated by comparing it to a pseudo-random binary 
sequence (PRBS) with the same amplitude constraints - a 
frequently-used input pattern for system identification A 
zero-mean Gaussian noise of variance 0.01 was added to y, and 
the output-error method in MATLAB System Identification 
Toolbox was applied. This test was performed for 500 times 
for each input signal with different noise patterns. Statistics 
of the resulting estimations are shown in Table [I] It is clear 
that substantially better parameter estimates are found with the 
proposed relaxation method. 
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the proposed 





0,1 


ai 


b 


True 


-1.8 


0.9 


0.1 


Mean (Opt) 


-1.800 


0.899 


0.10 


Std. Dev. (Opt) 


1.7xl0 -3 


1.7xl0" 3 


l.lxlO" 3 


Mean (PRBS) 


-1.499 


0.810 


0.0234 


Std. Dev. (PRBS) 


0.789 


0.383 


0.044 



TABLE I 

Parameter estimation comparison for the proposed 
optimization method and pseudo-random binary sequence 
(prbs) on example system {26}, 



VI. Conclusions 

We have proposed a new framework for design of input 
signals for system identification in the time domain. There are 
many potential applications in industrial processes, biomedical 
modelling, communication systems, etc. The time-domain 
input design problem is highly nonconvex. We propose using 
a convex relaxation based on semidefinite programming. 

For the case of a constraint on input amplitude, we have 
proven that the relaxation provides an upper bound which 
is accurate to within 2/tt. Furthermore we have provided 
a randomized strategy for finding a feasible solution which 
has a similar bound in terms of expected value of a reduced 
information matrix. For the case of a power constraint on the 
input, our method provides the true global optimum. It can 
also be extended in a natural way to constraints on the output, 
as well as multi-input multi-output and time-varying systems. 

A simple example illustrates the utility of the method, with 
the proposed strategy finding inputs which are far better than 
random binary sequences, resulting in substantial improve- 
ments in parameter estimates. 

As a direction for future work, it is well-known that a 
model which is optimal in terms of its parameter estimates 
may not be best for describing the response of the system 
under feedback control (see, e.g., 1251 . J4j) or long term open- 
loop simulation (see, e.g., ||26l , E71 ). It will be interesting to 
explore the application of the proposed method to problems 
of identification for control and simulation. 
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